#include "gbn_base_cuda.cuh"
using namespace std;

// 新增：确保实现文件中也能识别 Cairo（与头文件宏逻辑一致）
#ifdef GBN_ENABLE_CAIRO
#include <cairo/cairo.h>  // 再次包含，避免依赖顺序问题
// 必须包含 PDF 专用头文件（核心！很多人只包含 cairo.h 导致错误）
#include <cairo/cairo-pdf.h>  
#endif

// -------------------------- 数据结构定义 --------------------------
// -------------------------- Point点类型的实现 --------------------------
template <typename Float>
Point<Float>::Point(int dims) : coords(dims, 0.0f) {
    if (dims <= 0 || dims > MAX_DIMENSIONS) {
        throw std::invalid_argument(
            "Point: dimensions must be 1~" + std::to_string(MAX_DIMENSIONS) +
            " (current: " + std::to_string(dims) + ")"
        );
    }
}

template <typename Float>
Float& Point<Float>::operator[](int idx) {
    if (idx < 0 || idx >= (int)coords.size()) {
        throw std::out_of_range("Point: dimension index out of range");
    }
    return coords[idx];
}

template <typename Float>
const Float& Point<Float>::operator[](int idx) const {
    if (idx < 0 || idx >= (int)coords.size()) {
        throw std::out_of_range("Point: dimension index out of range");
    }
    return coords[idx];
}

template <typename Float>
int Point<Float>::getDims() const {
    return (int)coords.size();
}



// -------------------------- GPU核函数实现（高维平方距离计算） --------------------------
template <typename Float>
__global__ void computeDistSqKernel(int N, int dims, const Float* d_p, const Float* d_ref, Float* d_dist_sq) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= N) return;

    Float dist_sq = 0.0f;
    // 按维度循环累加平方距离（高维核心逻辑）
    for (int d = 0; d < dims; d++) {
        int idx = i * dims + d; // 扁平化数组索引：点i的第d维
        Float delta = d_p[idx] - d_ref[idx];
        dist_sq += delta * delta;
    }
    d_dist_sq[i] = dist_sq;
}

// -------------------------- 工具函数实现（随机数与GPU归约） --------------------------
template <typename Float>
Float randFloat(Float x) {
    if (x <= 0) {
        throw std::invalid_argument("randFloat: upper bound 'x' must be positive");
    }

    // 线程安全随机引擎（静态局部+线程局部存储）
    static thread_local std::mt19937_64 engine([]() {
        std::random_device rd;
        std::seed_seq seed{
            static_cast<unsigned int>(std::chrono::steady_clock::now().time_since_epoch().count()),
            rd(), rd(), rd(), rd()
        };
        return std::mt19937_64(seed);
    }());

    static thread_local std::uniform_real_distribution<Float> dist(0.0, 1.0);
    return dist(engine) * x;
}

// -------------------------- GPU核函数实现（高维平方距离计算） --------------------------
// 新增：GPU归约计算最大值的核函数实现（块内共享内存加速）
template <typename Float>
__global__ void gpuReduceMaxKernel(const Float* d_arr, int N, Float* d_temp) {
    // 1. 分配共享内存：存储当前线程块的局部元素（大小=线程块尺寸）
    __shared__ Float s_data[CUDA_BLOCK_SIZE];

    // 2. 计算线程全局索引，初始化共享内存
    const int tid = threadIdx.x;
    const int global_idx = blockIdx.x * blockDim.x + tid;

    // 初始化规则：覆盖数组元素的线程取原值，否则取类型最小值（不影响最大值计算）
    if (global_idx < N) {
        s_data[tid] = d_arr[global_idx];
    } else {
        // 用 std::numeric_limits 适配模板类型（比 FLT_MAX/DBL_MAX 更通用）
        s_data[tid] = std::numeric_limits<Float>::lowest();
    }

    // 3. 块内归约：按步长减半合并，计算块内最大值（共享内存加速）
    __syncthreads(); // 首次同步：确保所有线程完成共享内存初始化
    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            // 取当前线程与偏移线程的最大值（用 std::max 适配模板类型）
            s_data[tid] = std::max(s_data[tid], s_data[tid + s]);
        }
        __syncthreads(); // 每轮合并后同步，避免数据竞争
    }

    // 4. 块内结果写入临时缓冲区：仅块内第一个线程（tid=0）执行
    if (tid == 0) {
        d_temp[blockIdx.x] = s_data[0]; // 每个块的最大值存入 d_temp
    }
}

// 新增：主机端包装函数实现（调用核函数+错误处理+GPU-CPU同步）
template <typename Float>
Float gpuReduceMax(const Float* d_arr, int N, Float* d_temp) {
    // 1. 前置参数校验（避免无效计算或内存越界）
    if (N <= 0) {
        throw std::invalid_argument(
            "gpuReduceMax: array length N must be positive (current: " + std::to_string(N) + ")"
        );
    }
    if (d_arr == nullptr || d_temp == nullptr) {
        throw std::invalid_argument("gpuReduceMax: input GPU pointers (d_arr/d_temp) cannot be null");
    }

    // 2. 配置CUDA核函数参数（线程块=256，网格=向上取整覆盖所有元素）
    const dim3 block(CUDA_BLOCK_SIZE);
    const dim3 grid((N + block.x - 1) / block.x); // 网格大小=ceil(N / 256)

    // 检查临时缓冲区大小：d_temp 需至少容纳 grid.x 个元素（每个块的局部最大值）
    const size_t required_temp_bytes = grid.x * sizeof(Float);
    const size_t allocated_temp_bytes = CUDA_BLOCK_SIZE * sizeof(Float); // 类内 d_temp 分配大小
    if (allocated_temp_bytes < required_temp_bytes) {
        throw std::runtime_error(
            "gpuReduceMax: d_temp buffer size insufficient (required: " + std::to_string(required_temp_bytes) + 
            " bytes, allocated: " + std::to_string(allocated_temp_bytes) + " bytes)"
        );
    }

    // 3. 第一步归约：计算所有块的局部最大值，结果存入 d_temp（GPU端）
    gpuReduceMaxKernel<Float><<<grid, block>>>(d_arr, N, d_temp);
    // 检查核函数执行错误（复用类内的 checkCudaErr 逻辑思路）
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        throw std::runtime_error(
            "gpuReduceMax: first reduction kernel failed - " + std::string(cudaGetErrorString(err))
        );
    }
    // 同步GPU：确保第一步归约完成
    cudaDeviceSynchronize();
    if ((err = cudaGetLastError()) != cudaSuccess) {
        throw std::runtime_error(
            "gpuReduceMax: sync failed after first reduction - " + std::string(cudaGetErrorString(err))
        );
    }

    // 4. 第二步归约：若网格大小>1，对 d_temp 再次归约（直至只剩1个结果）
    Float h_result; // 最终结果（主机端）
    if (grid.x == 1) {
        // 优化：仅1个块时，直接读取 d_temp[0]
        cudaMemcpy(&h_result, d_temp, sizeof(Float), cudaMemcpyDeviceToHost);
    } else {
        // 二次归约：以 d_temp 为输入，网格大小=ceil(grid.x / 256)
        const dim3 grid2((grid.x + block.x - 1) / block.x);
        gpuReduceMaxKernel<Float><<<grid2, block>>>(d_temp, grid.x, d_temp);
        
        // 检查二次归约错误
        if ((err = cudaGetLastError()) != cudaSuccess) {
            throw std::runtime_error(
                "gpuReduceMax: second reduction kernel failed - " + std::string(cudaGetErrorString(err))
            );
        }
        cudaDeviceSynchronize();
        if ((err = cudaGetLastError()) != cudaSuccess) {
            throw std::runtime_error(
                "gpuReduceMax: sync failed after second reduction - " + std::string(cudaGetErrorString(err))
            );
        }

        // 读取最终结果
        cudaMemcpy(&h_result, d_temp, sizeof(Float), cudaMemcpyDeviceToHost);
    }

    // 5. 检查结果拷贝错误
    if ((err = cudaGetLastError()) != cudaSuccess) {
        throw std::runtime_error(
            "gpuReduceMax: copy result to CPU failed - " + std::string(cudaGetErrorString(err))
        );
    }

    return h_result;
}

// -------------------------- GBNBasePointSet 点集类成员函数实现 --------------------------
// 内部辅助函数：初始化
template <typename Float>
void GBNBasePointSet<Float>::init() {
    // 1. 初始化域范围（默认各维度为1.0）
    if (domain_ranges.empty()) {
        domain_ranges.resize(dimensions, 1.0f);
    } else if (domain_ranges.size() != (size_t)dimensions) {
        throw std::invalid_argument(
            "init: domain_ranges size (" + std::to_string(domain_ranges.size()) +
            ") must match dimensions (" + std::to_string(dimensions) + ")"
        );
    }
    for (auto& range : domain_ranges) {
        if (range <= 0) throw std::invalid_argument("init: domain range must be positive");
    }

    // 2. 计算归一化单位（统一各维度坐标到[0,1]，便于后续计算）
    Float max_domain = *std::max_element(domain_ranges.begin(), domain_ranges.end());
    unit = 1.0f / max_domain;

    // 3. 快照参数初始化（高维自动降级：禁用2D可视化格式）
    W = SNAP_VIS_SIZE;
    H = SNAP_VIS_SIZE;
    margin = 0.05f;
    accum_iterations = 0;
    recorded_frames = 0;
    frame_format = SNAP_FORMAT_TXT; // 高维默认仅保留TXT格式
    snap_cond = SnapshotCondition::POWER_OF_TWO; // 高维禁用「视觉变化触发」
    adjustFrameFormat(); // 自动过滤不支持的2D格式（PNG/PDF）

    // 4. 分配GPU内存（高维核心：扁平化数组存储，避免结构体内存对齐问题）
    size_t point_data_size = num_points * dimensions * sizeof(Float); // N×D 扁平化存储
    cudaMalloc(&d_points, point_data_size);
    cudaMalloc(&d_points_ref, point_data_size); // 参考点集（快照对比用）
    cudaMalloc(&d_dist_sq, num_points * sizeof(Float)); // 平方距离缓冲区
    cudaMalloc(&d_temp, CUDA_BLOCK_SIZE * sizeof(Float)); // 归约临时缓冲区
    checkCudaErr("init: GPU malloc failed");

    // 5. 初始化主机端点集（按各维度域范围随机分布）
    h_points.resize(num_points, Point<Float>(dimensions));
    for (auto& pt : h_points) {
        for (int d = 0; d < dimensions; d++) {
            // 各维度独立随机：[0, domain_ranges[d])
            pt[d] = randFloat(domain_ranges[d]);
        }
    }

    // 6. 同步主机端数据到GPU，并设置初始参考点集
    copyHostToDevice();
    copyDeviceToRef();

    // 7. 确保快照输出目录存在（若指定前缀）
    if (!snap_base_path.empty()) {
        std::string cmd = "mkdir -p " + snap_base_path;
        system(cmd.c_str());
    }
}

// 内部辅助函数：高维自动过滤2D格式（关键适配）
template <typename Float>
void GBNBasePointSet<Float>::adjustFrameFormat() {
    // 若维度≠2，强制过滤PNG/PDF格式，仅保留TXT
    if (dimensions != 2) {
        frame_format &= ~(SNAP_FORMAT_PNG  | SNAP_FORMAT_PDF);
        frame_format |= SNAP_FORMAT_TXT; // 确保TXT格式启用
        // 高维禁用「视觉变化触发」（无可视化，无法判断视觉变化）
        if (snap_cond == SnapshotCondition::VISUAL_CHANGE) {
            snap_cond = SnapshotCondition::POWER_OF_TWO;
            std::cout << "[Warning] High-dimensional mode: SnapshotCondition changed to POWER_OF_TWO (VISUAL_CHANGE disabled)" << std::endl;
        }
    }
}

// 内部辅助函数：检查CUDA错误（通用工具）
template <typename Float>
void GBNBasePointSet<Float>::checkCudaErr(const std::string& msg) const {
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        throw std::runtime_error(
            "CUDA Error [" + msg + "]: " + std::string(cudaGetErrorString(err))
        );
    }
}

// -------------------------- 构造函数实现 --------------------------
// 1. 随机生成高维点集（默认2D）
template <typename Float>
GBNBasePointSet<Float>::GBNBasePointSet(int num_points, 
                                        const std::vector<Float>& domain_ranges,
                                        const std::string& snap_path)
    : num_points(num_points), domain_ranges(domain_ranges), snap_base_path(snap_path) {
    // 合法性校验
    if (num_points <= 0) {
        throw std::invalid_argument("GBNBasePointSet: num_points must be positive (" + std::to_string(num_points) + ")");
    }
    if (num_points <= 0) {
        throw std::invalid_argument("GBNBasePointSet: num_points must be positive (" + std::to_string(num_points) + ")");
    }
    // 维度初始化：若域范围为空，默认2D；否则按域范围大小确定维度
    this->dimensions = domain_ranges.empty() ? 2 : static_cast<int>(domain_ranges.size());
    init(); // 调用内部初始化函数
}

// 2. 从TXT文件加载高维点集（格式：首行N D，后续N行每行D个坐标）
template <typename Float>
GBNBasePointSet<Float>::GBNBasePointSet(const std::string& txt_file_name, 
                                    const std::vector<Float>& domain_ranges,
                                    const std::string& snap_path)
                                    : domain_ranges(domain_ranges), snap_base_path(snap_path) {
    // 打开文件
    std::ifstream file(txt_file_name);
    if (!file.is_open()) {
        throw std::runtime_error("GBNBasePointSet: failed to open file - " + txt_file_name);
    }

    // 读取点数量与维度（TXT格式：首行 N D）
    int file_dims;
    file >> num_points >> file_dims;
    if (num_points <= 0 || file_dims <= 0) {
        throw std::invalid_argument("GBNBasePointSet: invalid N/D in file (N=" + std::to_string(num_points) + ", D=" + std::to_string(file_dims) + ")");
    }
    this->dimensions = file_dims;

    // 初始化GPU缓冲区
    init();

    // 读取高维坐标并赋值
    h_points.resize(num_points, Point<Float>(dimensions));
    for (int i = 0; i < num_points; i++) {
        if (file.eof()) {
            throw std::runtime_error("GBNBasePointSet: incomplete points in file (expected " + std::to_string(num_points) + ", got " + std::to_string(i) + ")");
        }
        for (int d = 0; d < dimensions; d++) {
            file >> h_points[i][d];
            // 校验坐标是否在域范围内
            if (h_points[i][d] < 0 || h_points[i][d] > domain_ranges[d]) {
                throw std::out_of_range("GBNBasePointSet: point " + std::to_string(i) + " dim " + std::to_string(d) + " out of range [0, " + std::to_string(domain_ranges[d]) + ")");
            }
        }
    }
    file.close();

    // 同步到GPU并更新参考点集
    copyHostToDevice();
    copyDeviceToRef();
}

// -------------------------- 域参数设置接口实现 --------------------------
template <typename Float>
void GBNBasePointSet<Float>::setToroidalMode(bool enable) {
    this->is_toroidal = enable;
    std::cout << "[Info] Toroidal mode " << (enable ? "enabled" : "disabled") << " (supports all dimensions)" << std::endl;
}

// -------------------------- 析构函数实现（释放GPU内存） --------------------------
template <typename Float>
GBNBasePointSet<Float>::~GBNBasePointSet() {
    // 安全释放GPU内存（避免野指针）
    if (d_points) cudaFree(d_points);
    if (d_points_ref) cudaFree(d_points_ref);
    if (d_dist_sq) cudaFree(d_dist_sq);
    if (d_temp) cudaFree(d_temp);
}

// -------------------------- 点集数据传输接口实现 --------------------------
// 1. 主机端 → GPU端（扁平化数组拷贝）
template <typename Float>
void GBNBasePointSet<Float>::copyHostToDevice() const{
    // 构建扁平化数组（主机端临时存储）
    std::vector<Float> h_flat(num_points * dimensions);
    for (int i = 0; i < num_points; i++) {
        for (int d = 0; d < dimensions; d++) {
            h_flat[i * dimensions + d] = h_points[i][d];
        }
    }
    // 拷贝到GPU
    cudaMemcpy(d_points, h_flat.data(), h_flat.size() * sizeof(Float), cudaMemcpyHostToDevice);
    checkCudaErr("copyHostToDevice: memcpy failed");
}

// 2. GPU端 → 主机端（还原为高维Point结构）
template <typename Float>
void GBNBasePointSet<Float>::copyDeviceToHost() const{
    // 从GPU拷贝扁平化数组
    std::vector<Float> h_flat(num_points * dimensions);
    cudaMemcpy(h_flat.data(), d_points, h_flat.size() * sizeof(Float), cudaMemcpyDeviceToHost);
    checkCudaErr("copyDeviceToHost: memcpy failed");
    // 还原为Point结构
    for (int i = 0; i < num_points; i++) {
        for (int d = 0; d < dimensions; d++) {
            h_points[i][d] = h_flat[i * dimensions + d];
        }
    }
}

// 3. GPU端 → GPU端（参考点集拷贝，避免CPU中转）
template <typename Float>
void GBNBasePointSet<Float>::copyDeviceToRef() const {
    size_t data_size = num_points * dimensions * sizeof(Float);
    cudaMemcpy(d_points_ref, d_points, data_size, cudaMemcpyDeviceToDevice);
    checkCudaErr("copyDeviceToRef: memcpy failed");
}

// -------------------------- 快照参数设置接口实现 --------------------------
template <typename Float>
void GBNBasePointSet<Float>::setOutBlackRatio(Float ratio) {
    if (ratio < 0 || ratio > 1) {
        throw std::invalid_argument("setOutBlackRatio: ratio must be [0,1] (" + std::to_string(ratio) + ")");
    }
    this->out_black_ratio = ratio;
    // 高维模式下警告（该参数仅2D有效）
    if (dimensions != 2) {
        std::cout << "[Warning] High-dimensional mode: out_black_ratio is only effective for 2D visualization" << std::endl;
    }
}

template <typename Float>
void GBNBasePointSet<Float>::setSnapshotCondition(SnapshotCondition cond) {
    // 高维禁用「视觉变化触发」
    if (dimensions != 2 && cond == SnapshotCondition::VISUAL_CHANGE) {
        std::cout << "[Warning] High-dimensional mode: VISUAL_CHANGE is disabled, using POWER_OF_TWO instead" << std::endl;
        snap_cond = SnapshotCondition::POWER_OF_TWO;
        return;
    }
    snap_cond = cond;
}

template <typename Float>
void GBNBasePointSet<Float>::setFrameBaseName(const std::string& base_name) {
    // 确保目录存在
    if (base_name.empty()) return;
    this->snap_base_path = base_name;
    
    // 自动创建地址对应的目录（提取路径中的目录部分，避免文件名干扰）
    size_t dir_pos = snap_base_path.find_last_of('/'); // Linux/macOS
    if (dir_pos == std::string::npos) {
        dir_pos = snap_base_path.find_last_of('\\'); // Windows
    }
    if (dir_pos != std::string::npos) {
        std::string dir_path = snap_base_path.substr(0, dir_pos);
        // 跨平台创建目录（忽略已存在的情况）
#ifdef _WIN32
        std::string cmd = "md \"" + dir_path + "\" 2>nul";
#else
        std::string cmd = "mkdir -p \"" + dir_path + "\"";
#endif
        system(cmd.c_str());
    }

    std::cout << "[Info] Snapshot base path updated to: " << snap_base_path << std::endl;
}

template <typename Float>
void GBNBasePointSet<Float>::setFrameFormat(int format_mask) {
    frame_format = format_mask;
    adjustFrameFormat(); // 应用高维格式过滤
}

template <typename Float>
void GBNBasePointSet<Float>::setFrameSize(int frame_w, int frame_h) {
    if (frame_w <= 0 || frame_h <= 0) {
        throw std::invalid_argument("setFrameSize: size must be positive (w=" + std::to_string(frame_w) + ", h=" + std::to_string(frame_h) + ")");
    }
    W = frame_w;
    H = frame_h;
    // 高维模式下警告（该参数仅2D有效）
    if (dimensions != 2) {
        std::cout << "[Warning] High-dimensional mode: frame size is only effective for 2D visualization" << std::endl;
    }
}

// -------------------------- 快照生成接口实现 --------------------------
/**
 * @brief 初始化快照默认参数（标准化配置，兼容高维/2D场景）
 * @note 1. 2D场景：启用TXT+PNG格式，按视觉变化触发快照
 *       2. 高维场景：自动降级为仅TXT格式，按2的幂次触发快照
 *       3. 统一默认输出路径、图像尺寸、边距等基础参数
 */
template <typename Float>
void GBNBasePointSet<Float>::initDefaultSnapshotParams() {
    // 若未通过构造函数设置地址（为空），自动赋值为SNAPSHOT_BASE_PATH
    if (snap_base_path.empty()) {
        this->setFrameBaseName(SNAPSHOT_BASE_PATH);
        std::cout << "[Info] Using default snapshot path: " << snap_base_path << std::endl;
    }

    // 2. 默认快照格式（根据维度自动适配：2D多格式，高维仅TXT）
    int default_format = 0;
    if (this->dimensions == 2) {
        // 2D场景：启用TXT（数据存储）+ PNG（快速可视化）
        default_format = SNAP_FORMAT_TXT | SNAP_FORMAT_PNG;
    } else {
        // 高维场景：仅保留TXT格式（高维无可视化需求，TXT便于后续数据加载/分析）
        default_format = SNAP_FORMAT_TXT;
        std::cout << "[Info] High-dimensional mode (" << this->dimensions << "D): Default snapshot format set to TXT only" << std::endl;
    }
    this->setFrameFormat(default_format); // 调用格式设置接口（自动过滤不支持的格式）

    // 3. 默认快照触发条件（根据维度适配：2D视觉变化，高维幂次迭代）
    SnapshotCondition default_cond = SnapshotCondition::VISUAL_CHANGE;
    if (this->dimensions != 2) {
        // 高维禁用"视觉变化"（无图像像素概念），改用"2的幂次"（迭代1,2,4,8...时生成，平衡数据量与趋势观察）
        default_cond = SnapshotCondition::POWER_OF_TWO;
    }
    this->setSnapshotCondition(default_cond); // 调用条件设置接口（自动处理高维兼容）

    // 4. 2D默认可视化参数（高维调用时会自动警告，无需额外判断）
    this->setFrameSize(800, 600); // 默认图像尺寸：800x600（兼顾清晰度与存储大小）
    this->margin = 0.03f;         // 默认边距：3%（避免点贴边，提升可视化效果）
    this->setOutBlackRatio(0.02f); // 默认黑度比例：2%（控制点大小，避免重叠/过小）

    // 5. 初始化快照计数（确保从0开始，避免统计错误）
    this->accum_iterations = 0;
    this->recorded_frames = 0;

    std::cout << "[Info] Default snapshot parameters initialized successfully:" << std::endl;
    std::cout << "  - Output Path: " << this->snap_base_path << "frame_XXXXXXX.{txt/png/PDF}" << std::endl;
    std::cout << "  - Trigger Condition: " << (default_cond == SnapshotCondition::VISUAL_CHANGE ? "Visual Change (2D)" : "Power of Two (High-D)") << std::endl;
    std::cout << "  - Format: " << (default_format == SNAP_FORMAT_TXT ? "TXT" : "TXT + PNG") << std::endl;
}

template <typename Float>
void GBNBasePointSet<Float>::generateSnapshot(int iter, bool force) {
    accum_iterations++;
    bool need_snap = force;

    // 按触发条件判断是否生成快照
    if (!need_snap) {
        switch (snap_cond) {
            case SnapshotCondition::VISUAL_CHANGE: {
                // 仅2D有效：计算点移动是否超过1像素
                Float max_pixel_dist = calcMaxPixelDist();
                if (max_pixel_dist >= 1.0f) {
                    need_snap = true;
                    copyDeviceToRef(); // 更新参考点集
                    std::cout << "[Snapshot] Iter " << iter << " (visual change: " << max_pixel_dist << "px)" << std::endl;
                }
                break;
            }
            case SnapshotCondition::POWER_OF_TWO: {
                // 迭代次数为2的幂（1,2,4,8...）
                need_snap = (iter > 0) && ((iter & (iter - 1)) == 0);
                break;
            }
            case SnapshotCondition::EVERY_ITER: {
                // 每次迭代都生成
                need_snap = true;
                break;
            }
        }
    }

    // 无需快照或无有效格式，直接返回
    if (!need_snap || frame_format == 0) return;

    // 生成带前导零的文件名前缀（如 "frame_0000100"）
    char iter_str[8];
    sprintf(iter_str, "%07d", iter);
    std::string frame_prefix = "frame_" + std::string(iter_str);
    recorded_frames++;

    // 同步GPU数据到主机端（用于保存）
    copyDeviceToHost();

    // -------------------------- 关键修改：格式-文件夹映射与路径生成 --------------------------
    // 3.1 定义格式对应的子文件夹名称（可根据需求调整）
    const std::unordered_map<int, std::string> format_folder_map = {
        {SNAP_FORMAT_TXT,  "frametxt"},  // TXT格式 → 子文件夹 frametxt
        {SNAP_FORMAT_PNG,  "framepng"},  // PNG格式 → 子文件夹 framepng
        {SNAP_FORMAT_PDF,  "framepdf"}   // PDF格式 → 子文件夹 framepdf
    };

    // 3.2 遍历所有启用的格式，生成对应路径并保存
    for (const auto& [format, folder] : format_folder_map) {
        // 跳过当前未启用的格式
        if (!(frame_format & format)) continue;

        // 3.2.1 生成格式专属的完整路径（基础路径 + 子文件夹 + 文件名）
        std::string format_base_path = snap_base_path;
        // 确保基础路径末尾有路径分隔符（避免 "outputframetxt" 这类错误）
        if (!format_base_path.empty() && format_base_path.back() != '/' && format_base_path.back() != '\\') {
#ifdef _WIN32
            format_base_path += "\\"; // Windows路径分隔符
#else
            format_base_path += "/";  // Linux/macOS路径分隔符
#endif
        }
        // 拼接子文件夹（如 "output/" + "frametxt/"）
        const std::string format_folder_path = format_base_path + folder + "/";
        // 拼接完整文件路径（如 "output/frametxt/frame_0000100.txt"）
        const std::string file_ext = (format == SNAP_FORMAT_TXT) ? ".txt" : 
                                     (format == SNAP_FORMAT_PNG) ? ".png" : ".pdf";
        const std::string full_file_path = format_folder_path + frame_prefix + file_ext;

        // 3.2.2 自动创建子文件夹（跨平台兼容，忽略已存在的情况）
#ifdef _WIN32
        // Windows：md 命令创建目录，2>nul 忽略"目录已存在"错误
        const std::string mkdir_cmd = "md \"" + format_folder_path + "\" 2>nul";
#else
        // Linux/macOS：mkdir -p 递归创建目录（无则创建，有则忽略）
        const std::string mkdir_cmd = "mkdir -p \"" + format_folder_path + "\"";
#endif
        system(mkdir_cmd.c_str()); // 执行创建目录命令

        // 3.2.3 按格式保存文件（调用原保存接口，路径改为格式专属路径）
        try {
            if (format == SNAP_FORMAT_TXT) {
                saveTXT(full_file_path);
            } else if (format == SNAP_FORMAT_PNG) {
                savePNG(full_file_path);
            } else if (format == SNAP_FORMAT_PDF) {
                savePDF(full_file_path);
            }
            std::cout << "[Snapshot Saved] " << full_file_path << std::endl; // 打印保存路径，便于调试
        } catch (const std::exception& e) {
            // 捕获单格式保存失败的异常，避免影响其他格式（可选增强）
            std::cerr << "[Warning] Failed to save " << file_ext << " snapshot: " << e.what() << std::endl;
        }
    }
/*
    // 按格式生成文件
    if (frame_format & SNAP_FORMAT_TXT) {
        saveTXT(frame_path + ".txt");
    }
    if (frame_format & SNAP_FORMAT_PNG) {
        savePNG(frame_path + ".png");
    }
    if (frame_format & SNAP_FORMAT_PDF) {
        savePDF(frame_path + ".pdf");
    }
*/
    accum_iterations = 0;
}

// -------------------------- 多格式保存接口实现 --------------------------
/**
 * @brief 保存高维点集为TXT文件（支持后续加载，兼容任意维度）
 * @param file_path 输出文件路径（如"output/frame_0000100.txt"）
 * @throw std::runtime_error 当文件打开失败时抛出异常
 * @note TXT格式规范：
 *       - 首行：2个整数，分别为「点集数量 N」和「维度 D」
 *       - 第2~N+1行：每行 D 个浮点数，对应1个点的各维度坐标（已归一化到[0,1]）
 *       - 坐标精度：6位小数，确保数据精度与加载一致性
 */
template <typename Float>
void GBNBasePointSet<Float>::saveTXT(const std::string& file_path) const {
    // 打开文件（写模式，若文件不存在则创建，存在则覆盖）
    std::ofstream file(file_path);
    if (!file.is_open()) {
        throw std::runtime_error(
            "saveTXT: failed to open file - " + file_path + 
            " (check path validity or file permissions)"
        );
    }

    // 1. 写入文件头：点数量 + 维度（高维核心标识）
    file << num_points << " " << dimensions << std::endl;

    // 2. 配置输出精度（6位小数，平衡精度与文件大小）
    file << std::fixed << std::setprecision(6);

    // 3. 遍历所有点，写入各维度归一化坐标（统一映射到[0,1]域）
    for (const auto& pt : h_points) {
        // 遍历当前点的所有维度（从0到dimensions-1）
        for (int d = 0; d < dimensions; d++) {
            // 坐标归一化：原始坐标 × 单位系数（unit = 1/最大域范围）
            Float normalized_coord = pt[d] * unit;
            // 维度间用空格分隔，最后一个维度后换行
            file << normalized_coord << (d == dimensions - 1 ? "\n" : " ");
        }
    }

    // 关闭文件（确保缓冲区数据刷写完成）
    file.close();
}

/**
 * @brief 保存2D点集为PNG图像（仅2D有效，高维调用时会触发警告）
 * @param file_path 输出文件路径（如"output/frame_0000100.png"）
 * @throw std::runtime_error 当文件保存失败或非2D模式时抛出异常
 */
template <typename Float>
void GBNBasePointSet<Float>::savePNG(const std::string& file_path) const {
    // 非2D模式直接报错（PNG仅支持视觉化2D点集）
    if (dimensions != 2) {
        throw std::runtime_error(
            "savePNG: only 2D point sets are supported (current dimensions: " + 
            std::to_string(dimensions) + ")"
        );
    }

    // 1. 创建白色背景图像（OpenCV默认通道顺序：BGR）
    cv::Mat img(H, W, CV_8UC3, SNAP_BG_COLOR);
    if (img.empty()) {
        throw std::runtime_error(
            "savePNG: failed to create image - invalid size (W=" + 
            std::to_string(W) + ", H=" + std::to_string(H) + ")"
        );
    }

    // 2. 计算点绘制半径（适配图像尺寸，避免点重叠或过小）
    // 逻辑：基于默认尺寸(SNAP_VIS_SIZE)的比例缩放，最小半径为1像素
    Float point_rad = SNAP_POINT_RAD<Float> * std::min(W, H) / SNAP_VIS_SIZE;
    point_rad = std::max(static_cast<Float>(1.0), point_rad); // 确保不小于1像素

    // 3. 遍历2D点集，绘制每个点（域坐标→图像坐标转换）
    for (const auto& pt : h_points) {
        // 提取2D点的x/y坐标（高维已提前过滤，此处直接取前2维）
        Float x_domain = pt[0] * unit; // 归一化到[0,1]的x域坐标
        Float y_domain = pt[1] * unit; // 归一化到[0,1]的y域坐标

        // 坐标映射：域坐标 → 图像像素坐标（处理OpenCV左上角原点与边距）
        // 公式：像素坐标 = (域坐标 + 边距) × 图像尺寸 / (1 + 2×边距)
        // 注：y轴翻转（因OpenCV图像原点在左上角，域原点在左下角）
        int x_pixel = static_cast<int>(
            (x_domain + margin) * W / (1 + 2 * margin)
        );
        int y_pixel = static_cast<int>(
            (1 - (y_domain + margin)) * H / (1 + 2 * margin)
        );

        // 确保像素坐标在图像范围内（避免超出图像边界导致绘制失败）
        x_pixel = std::clamp(x_pixel, 0, W - 1);
        y_pixel = std::clamp(y_pixel, 0, H - 1);

        // 绘制红色填充圆（参数：图像、像素坐标、半径、颜色、填充方式）
        cv::circle(
            img, 
            cv::Point(x_pixel, y_pixel), 
            static_cast<int>(point_rad), 
            SNAP_POINT_COLOR, 
            -1 // 填充圆（-1=填充，≥0=仅绘制边框）
        );
    }

    // 4. 保存PNG文件（压缩等级3：平衡保存速度与文件体积）
    std::vector<int> png_params = {cv::IMWRITE_PNG_COMPRESSION, 3};
    if (!cv::imwrite(file_path, img, png_params)) {
        throw std::runtime_error(
            "savePNG: failed to save image - " + file_path + 
            " (check OpenCV PNG encoder or disk space)"
        );
    }
}

/**
 * @brief 保存2D点集为PDF矢量文件（仅2D有效，支持打印与多页）
 * @param file_path 输出文件路径（如"output/frame_0000100.pdf"）
 * @throw std::runtime_error 当PDF创建失败或非2D模式时抛出异常
 */
template <typename Float>
void GBNBasePointSet<Float>::savePDF(const std::string& file_path) const {
    // 非2D模式直接报错（PDF可视化仅支持2D）
    if (dimensions != 2) {
        throw std::runtime_error(
            "savePDF: only 2D point sets are supported (current dimensions: " + 
            std::to_string(dimensions) + ")"
        );
    }

    // 1. 创建PDF表面（A4纸张尺寸：595pt × 842pt，符合PDF标准）
    cairo_surface_t* pdf_surface = cairo_pdf_surface_create(
        file_path.c_str(),  // 输出路径
        595.0,              // A4宽度（pt，1pt = 1/72英寸）
        842.0               // A4高度（pt）
    );
    if (cairo_surface_status(pdf_surface) != CAIRO_STATUS_SUCCESS) {
        throw std::runtime_error(
            "savePDF: failed to create PDF surface - " + file_path + 
            " (check cairo PDF backend or disk space)"
        );
    }

    // 2. 创建Cairo绘图上下文（绑定PDF表面）
    cairo_t* cr = cairo_create(pdf_surface);
    if (cairo_status(cr) != CAIRO_STATUS_SUCCESS) {
        cairo_surface_destroy(pdf_surface); // 释放已创建的表面
        throw std::runtime_error("savePDF: failed to create cairo context");
    }

    // 3. 配置PDF绘图参数（坐标变换、颜色、点半径）
    // 3.1 坐标变换：应用A4页边距（50pt），适配点集域范围
    const Float page_margin = 50.0f; // A4页边距（避免内容超出打印区域）
    const Float content_width = 595.0f - 2 * page_margin;  // 内容区宽度
    const Float content_height = 842.0f - 2 * page_margin; // 内容区高度

    // 平移到内容区起点 + 缩放适配内容区（统一映射域坐标到PDF内容区）
    cairo_translate(cr, page_margin, page_margin);
    cairo_scale(
        cr, 
        content_width / (static_cast<Float>(W) * (1 + 2 * margin)),  // x轴缩放
        content_height / (static_cast<Float>(H) * (1 + 2 * margin))  // y轴缩放
    );
    cairo_translate(cr, static_cast<Float>(W) * margin, static_cast<Float>(H) * margin); // 应用点集边距

    // 3.2 设置点颜色（黑色填充，RGB: 0,0,0）
    cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);

    // 3.3 计算点半径（基于黑度比例）
    Float total_black_area = static_cast<Float>(W * H) * this->out_black_ratio;
    Float point_rad = sqrt(total_black_area / (num_points * PI<Float>));

    // 4. 遍历2D点集，绘制PDF矢量点
    for (const auto& pt : h_points) {
        // 计算点的PDF坐标（归一化后适配缩放后的内容区）
        Float x_pdf = pt[0] * unit * static_cast<Float>(W);
        Float y_pdf = pt[1] * unit * static_cast<Float>(H);

        // 绘制填充圆（Cairo指令：圆弧→闭合路径→填充）
        cairo_arc(
            cr, 
            x_pdf,      // 圆心x坐标
            y_pdf,      // 圆心y坐标
            point_rad,  // 半径
            0.0,        // 起始角度（弧度，0=右向）
            2 * PI<Float> // 结束角度（完整圆）
        );
        cairo_fill(cr); // 填充圆（使用当前颜色）
    }

    // 5. 释放Cairo资源（避免内存泄漏）
    cairo_show_page(cr);          // 刷新PDF页面（确保内容写入）
    cairo_destroy(cr);            // 销毁绘图上下文
    cairo_surface_destroy(pdf_surface); // 销毁PDF表面（触发文件写入）
}

// -------------------------- 辅助计算接口实现（高维兼容+模板化） --------------------------
/**
 * @brief 计算点集与参考点集的最大像素距离（仅2D有效，用于视觉变化判断）
 * @return 最大像素距离（2D返回实际值，高维返回0并打印警告）
 * @throw std::runtime_error 当CUDA核函数执行失败时抛出异常
 */
template <typename Float>
Float GBNBasePointSet<Float>::calcMaxPixelDist() {
    // 高维模式：直接返回0并警告（无可视化像素概念）
    if (dimensions != 2) {
        std::cout << "[Warning] High-dimensional mode: calcMaxPixelDist is only effective for 2D (return 0)" << std::endl;
        return 0.0f;
    }

    // 1. 配置CUDA核函数执行参数（线程块/网格大小）
    dim3 block(CUDA_BLOCK_SIZE);
    dim3 grid((num_points + block.x - 1) / block.x); // 向上取整确保覆盖所有点

    // 2. 调用GPU核函数计算每个点的平方距离（使用模板化核函数computeDistSqKernel）
    computeDistSqKernel<Float><<<grid, block>>>(
        num_points,   // 点集数量
        dimensions,   // 维度（固定2D）
        d_points,     // 当前点集（GPU）
        d_points_ref, // 参考点集（GPU，用于对比变化）
        d_dist_sq     // 输出：平方距离数组（GPU）
    );
    checkCudaErr("calcMaxPixelDist: computeDistSqKernel execution failed");

    // 3. GPU归约计算最大平方距离，转换为实际距离
    Float max_dist_sq = gpuReduceMax<Float>(d_dist_sq, num_points, d_temp);
    Float max_dist = sqrt(max_dist_sq);

    // 4. 将域距离转换为像素距离（单位域尺寸 → 像素数映射）
    // 公式：像素距离 = 域距离 × 单位系数 × 有效像素密度
    Float effective_pixel_density = std::min(W, H) / (1 + 2 * margin); // 去除边距后的像素密度
    Float max_pixel_dist = max_dist * unit * effective_pixel_density;

    return max_pixel_dist;
}

/**
 * @brief 计算点集平均最近邻距离（支持高维，评估点集分布均匀性）
 * @return 归一化后的平均最近邻距离（范围[0,1]）
 * @note 高维场景下自动适配维度，2D场景兼容原逻辑；小数据量用CPU实现，大数据量可扩展GPU加速
 */
template <typename Float>
Float GBNBasePointSet<Float>::calcAverageDist() const {
    // 1. 若点集数量≤1，无最近邻可计算（返回0避免除以0）
    if (num_points <= 1) {
        std::cout << "[Warning] calcAverageDist: num_points ≤ 1 (return 0)" << std::endl;
        return 0.0f;
    }

    // 2. 将GPU端数据同步到主机端（CPU计算需访问点坐标）
    copyDeviceToHost();

    Float total_min_dist = 0.0f; // 所有点的最近邻距离总和
    Float max_possible_dist = 0.0f; // 最大可能距离（用于归一化，高维适配）

    // 3. 计算最大可能距离（高维欧几里得距离，基于各维度域范围）
    for (int d = 0; d < dimensions; d++) {
        Float dim_range = domain_ranges[d] * unit; // 归一化后的维度范围
        max_possible_dist += dim_range * dim_range;
    }
    max_possible_dist = sqrt(max_possible_dist); // 高维空间对角线距离

    // 4. 遍历每个点，寻找其最近邻（高维适配：遍历所有维度计算距离）
    for (int i = 0; i < num_points; i++) {
        Float min_dist = max_possible_dist; // 初始化为最大可能距离（确保能找到更小值）

        // 对比当前点与其他所有点（j≠i）
        for (int j = 0; j < num_points; j++) {
            if (i == j) continue; // 跳过自身

            // 4.1 计算高维欧几里得距离（累加各维度平方差）
            Float dist_sq = 0.0f;
            for (int d = 0; d < dimensions; d++) {
                Float delta = h_points[i][d] * unit - h_points[j][d] * unit; // 归一化后计算差值
                dist_sq += delta * delta;
            }
            Float dist = sqrt(dist_sq);

            // 4.2 更新当前点的最近邻距离
            if (dist < min_dist) {
                min_dist = dist;
            }
        }

        // 4.3 累加当前点的最近邻距离（处理环形域场景：若启用环形，需额外计算边界距离）
        if (this->is_toroidal) {
            // 环形域特殊处理：计算点到"边界对端"的距离（如2D中x=0与x=W视为相邻）
            Float toroidal_min_dist = min_dist;
            for (int j = 0; j < num_points; j++) {
                if (i == j) continue;

                Float toroidal_dist_sq = 0.0f;
                for (int d = 0; d < dimensions; d++) {
                    Float dim_range = domain_ranges[d] * unit;
                    Float delta = fabs(h_points[i][d] * unit - h_points[j][d] * unit);
                    // 环形域中，距离取"直接距离"与"绕边界距离"的最小值
                    delta = std::min(delta, dim_range - delta);
                    toroidal_dist_sq += delta * delta;
                }
                Float toroidal_dist = sqrt(toroidal_dist_sq);
                if (toroidal_dist < toroidal_min_dist) {
                    toroidal_min_dist = toroidal_dist;
                }
            }
            total_min_dist += toroidal_min_dist;
        } else {
            total_min_dist += min_dist;
        }
    }

    // 5. 计算平均最近邻距离（归一化到[0,1]范围，便于跨场景对比）
    Float average_min_dist = total_min_dist / num_points;
    Float normalized_average_dist = average_min_dist / max_possible_dist;

    return normalized_average_dist;
}

/**
 * 计算点集的径向功率谱（支持高维，评估点集分布的频率特性）
 * @param[out] spectrum 输出：径向功率谱数组（索引对应径向距离区间，值为该区间的功率密度）
 * @param[in] num_bins 径向距离区间数量（默认100，区间数越多分辨率越高）
 * @param[in] max_radius 最大径向距离（默认0.0，自动计算为高维空间对角线的1/2）
 * @throw invalid_argument 当 num_bins ≤0 或 max_radius <0 时抛出
 * @note 1. 高维场景：基于点对欧氏距离统计径向分布，功率定义为"单位径向区间内的点对数量占比"；
 *       2. 环形域场景：自动修正距离计算（取最短环形距离）；
 *       3. 结果已归一化：所有区间功率之和为1，便于跨点集对比。
 */
template <typename Float>
void GBNBasePointSet<Float>::calcRadialPowerSpectrum(std::vector<Float>& spectrum, int num_bins, Float max_radius) {
    // 1. 参数合法性校验
    if (num_bins <= 0) {
        throw std::invalid_argument(
            "calcRadialPowerSpectrum: num_bins must be positive (current: " + std::to_string(num_bins) + ")"
        );
    }
    if (max_radius < 0.0f) {
        throw std::invalid_argument(
            "calcRadialPowerSpectrum: max_radius cannot be negative (current: " + std::to_string(max_radius) + ")"
        );
    }

    // 2. 初始化输出功率谱数组（清零）
    spectrum.resize(num_bins, 0.0f);

    // 3. 计算总点对数（避免重复统计：i<j，共 N*(N-1)/2 对）
    const int total_pairs = num_points * (num_points - 1) / 2;
    if (total_pairs <= 0) {
        std::cout << "[Warning] calcRadialPowerSpectrum: too few points (num_points ≤1), spectrum remains zero" << std::endl;
        return;
    }

    // 4. 自动计算最大径向距离（若未指定）
    // 高维空间最大距离：归一化后各维度范围均为[0,1]，对角线长度=sqrt(dimensions)，取1/2作为最大径向半径
    if (max_radius <= 1e-6f) {
        max_radius = static_cast<Float>(sqrt(dimensions)) / 2.0f;
    }
    const Float bin_size = max_radius / num_bins; // 每个径向区间的宽度

    // 5. 遍历所有点对，统计径向距离分布（核心逻辑）
    for (int i = 0; i < num_points; i++) {
        for (int j = i + 1; j < num_points; j++) {
            // 5.1 计算点对的欧氏距离（高维适配，支持环形域修正）
            Float dist_sq = 0.0f;
            for (int d = 0; d < dimensions; d++) {
                Float delta = h_points[i][d] * unit - h_points[j][d] * unit; // 归一化到[0,1]域
                
                // 环形域修正：取该维度的最短距离（如x=0.9与x=0.1的距离=0.2，而非0.8）
                if (is_toroidal) {
                    const Float dim_range = 1.0f; // 归一化后维度范围为[0,1]
                    delta = fabs(delta);
                    delta = std::min(delta, dim_range - delta);
                }
                
                dist_sq += delta * delta;
            }
            const Float dist = sqrt(dist_sq);

            // 5.2 跳过超出最大径向半径的点对
            if (dist >= max_radius) {
                continue;
            }

            // 5.3 确定当前距离所属的径向区间，累加计数
            const int bin_idx = static_cast<int>(dist / bin_size);
            if (bin_idx >= 0 && bin_idx < num_bins) {
                spectrum[bin_idx] += 1.0f;
            }
        }
    }

    // 6. 功率谱归一化：每个区间的功率 = 区间点对数 / 总点对数（确保总和为1）
    for (auto& power : spectrum) {
        power /= total_pairs;
    }

    std::cout << "[Info] Radial power spectrum computed: " << num_bins << " bins, max radius = " << max_radius << std::endl;
}

/**
 * 显示径向功率谱（2D 可视化，仅支持 2D 场景，高维场景触发警告）
 * @param[in] spectrum 输入：已计算的径向功率谱数组
 * @param[in] window_name 显示窗口名称（默认"Radial Power Spectrum"）
 * @param[in] plot_size 显示图像尺寸（默认512x512，仅2D有效）
 * @throw runtime_error 当非2D场景或功率谱为空时抛出
 * @note 1. 可视化内容：x轴为"归一化径向距离"，y轴为"功率密度"，红色曲线绘制功率谱；
 *       2. 支持窗口交互（按任意键关闭窗口）；
 *       3. 自动添加坐标轴标签与网格，提升可读性。
 */
 template <typename Float>
void GBNBasePointSet<Float>::showRadialPowerSpectrum(const std::vector<Float>& spectrum, 
                                                     const std::string& window_name, 
                                                     int plot_size) {
    // 1. 合法性校验
    if (dimensions != 2) {
        throw std::runtime_error(
            "showRadialPowerSpectrum: only 2D point sets are supported (current dimensions: " + 
            std::to_string(dimensions) + ")"
        );
    }
    if (spectrum.empty()) {
        throw std::runtime_error("showRadialPowerSpectrum: input spectrum is empty (call calcRadialPowerSpectrum first)");
    }
    if (plot_size <= 0) {
        throw std::invalid_argument(
            "showRadialPowerSpectrum: plot_size must be positive (current: " + std::to_string(plot_size) + ")"
        );
    }

    // 2. 创建白色背景图像（用于绘制功率谱）
    const int margin = 50; // 图像边距（预留坐标轴标签空间）
    const int plot_w = plot_size;
    const int plot_h = plot_size;
    cv::Mat img(plot_h, plot_w, CV_8UC3, SNAP_BG_COLOR); // 白色背景

    // 3. 计算绘图参数（适配图像尺寸）
    const int num_bins = static_cast<int>(spectrum.size());
    const Float max_power = *std::max_element(spectrum.begin(), spectrum.end()); // 最大功率（用于y轴归一化）
    const Float bin_width = static_cast<Float>(plot_w - 2 * margin) / num_bins; // 每个区间的绘图宽度
    const Float y_scale = static_cast<Float>(plot_h - 2 * margin) / (max_power > 1e-6f ? max_power : 1.0f); // y轴缩放系数

    // 4. 绘制功率谱曲线（红色，线宽2）
    for (int bin_idx = 0; bin_idx < num_bins - 1; bin_idx++) {
        // 计算当前区间与下一个区间的绘图坐标
        const int x1 = margin + static_cast<int>(bin_idx * bin_width);
        const int x2 = margin + static_cast<int>((bin_idx + 1) * bin_width);
        // y轴翻转（OpenCV图像原点在左上角，需将"功率高"的点画在上方）
        const int y1 = plot_h - margin - static_cast<int>(spectrum[bin_idx] * y_scale);
        const int y2 = plot_h - margin - static_cast<int>(spectrum[bin_idx + 1] * y_scale);

        // 绘制线段（红色，线宽2）
        cv::line(img, cv::Point(x1, y1), cv::Point(x2, y2), cv::Scalar(0, 0, 255), 2);
    }

    // 5. 添加坐标轴标签与网格
    // 5.1 坐标轴标签（黑色，字体大小0.5）
    const cv::Scalar text_color(0, 0, 0);
    cv::putText(img, "Normalized Radial Distance", cv::Point(margin, plot_h - 10), 
                cv::FONT_HERSHEY_SIMPLEX, 0.5, text_color);
    cv::putText(img, "Power Density", cv::Point(10, margin), 
                cv::FONT_HERSHEY_SIMPLEX, 0.5, text_color, 1, cv::LINE_AA, true); // 垂直文字

    // 5.2 网格线（灰色，线宽1，虚线）
    const cv::Scalar grid_color(200, 200, 200);
    const int grid_step = 5; // 每5个区间画一条竖线
    for (int bin_idx = 0; bin_idx < num_bins; bin_idx += grid_step) {
        const int x = margin + static_cast<int>(bin_idx * bin_width);
        cv::line(img, cv::Point(x, margin), cv::Point(x, plot_h - margin), grid_color, 1, cv::LINE_AA);
    }

    // 6. 显示窗口并等待交互
    cv::namedWindow(window_name, cv::WINDOW_NORMAL);
    cv::imshow(window_name, img);
    cv::resizeWindow(window_name, plot_w, plot_h);
    std::cout << "[Info] Press any key to close " << window_name << " window..." << std::endl;
    cv::waitKey(0); // 等待按键关闭窗口
    cv::destroyWindow(window_name);
}

/**
 * 保存径向功率谱图像（2D 矢量图/位图，仅支持2D场景）
 * @param[in] spectrum 输入：已计算的径向功率谱数组
 * @param[in] file_path 输出文件路径（支持.png/.pdf格式）
 * @param[in] plot_size 图像尺寸（默认800x600）
 * @throw runtime_error 当非2D场景、功率谱为空或文件保存失败时抛出
 */
template <typename Float>
void GBNBasePointSet<Float>::saveRadialPowerSpectrum(const std::vector<Float>& spectrum, 
                                                     const std::string& file_path, 
                                                     int plot_size) {
    // 1. 合法性校验（同显示函数）
    if (dimensions != 2) {
        throw std::runtime_error(
            "saveRadialPowerSpectrum: only 2D point sets are supported (current dimensions: " + 
            std::to_string(dimensions) + ")"
        );
    }
    if (spectrum.empty()) {
        throw std::runtime_error("saveRadialPowerSpectrum: input spectrum is empty (call calcRadialPowerSpectrum first)");
    }
    if (plot_size <= 0) {
        throw std::invalid_argument(
            "saveRadialPowerSpectrum: plot_size must be positive (current: " + std::to_string(plot_size) + ")"
        );
    }

    // 2. 绘制功率谱图像（复用显示函数的绘图逻辑）
    const int margin = 50;
    const int plot_w = plot_size;
    const int plot_h = plot_size;
    cv::Mat img(plot_h, plot_w, CV_8UC3, SNAP_BG_COLOR);

    const int num_bins = static_cast<int>(spectrum.size());
    const Float max_power = *std::max_element(spectrum.begin(), spectrum.end());
    const Float bin_width = static_cast<Float>(plot_w - 2 * margin) / num_bins;
    const Float y_scale = static_cast<Float>(plot_h - 2 * margin) / (max_power > 1e-6f ? max_power : 1.0f);

    // 绘制功率谱曲线
    for (int bin_idx = 0; bin_idx < num_bins - 1; bin_idx++) {
        const int x1 = margin + static_cast<int>(bin_idx * bin_width);
        const int x2 = margin + static_cast<int>((bin_idx + 1) * bin_width);
        const int y1 = plot_h - margin - static_cast<int>(spectrum[bin_idx] * y_scale);
        const int y2 = plot_h - margin - static_cast<int>(spectrum[bin_idx + 1] * y_scale);
        cv::line(img, cv::Point(x1, y1), cv::Point(x2, y2), cv::Scalar(0, 0, 255), 2);
    }

    // 添加坐标轴标签
    const cv::Scalar text_color(0, 0, 0);
    cv::putText(img, "Normalized Radial Distance", cv::Point(margin, plot_h - 10), 
                cv::FONT_HERSHEY_SIMPLEX, 0.5, text_color);
    cv::putText(img, "Power Density", cv::Point(10, margin), 
                cv::FONT_HERSHEY_SIMPLEX, 0.5, text_color, 1, cv::LINE_AA, true);

    // 3. 根据文件后缀保存对应格式
    const std::string ext = file_path.substr(file_path.find_last_of('.') + 1);
    if (ext == "png") {
        std::vector<int> png_params = {cv::IMWRITE_PNG_COMPRESSION, 3};
        if (!cv::imwrite(file_path, img, png_params)) {
            throw std::runtime_error("saveRadialPowerSpectrum: failed to save PNG - " + file_path);
        }
    } else if ( ext == "pdf") {
        // 复用原有 savePDF 的 Cairo 绘图逻辑（简化版）
        cairo_surface_t* surface = cairo_pdf_surface_create(file_path.c_str(), plot_w, plot_h);
        if (!surface) {
            throw std::runtime_error("saveRadialPowerSpectrum: failed to create " + ext + " surface - " + file_path);
        }
    
        // 创建 Cairo 绘图上下文（绑定 PDF 表面）
        cairo_t* cr = cairo_create(surface);
        if (cairo_status(cr) != CAIRO_STATUS_SUCCESS) {
        // 1. 先销毁表面（无论格式，避免内存泄漏）
        cairo_surface_destroy(surface);
        
        // 2. 仅 EPS 需要手动关闭文件句柄，但当前已移除 EPS 支持，直接删除该逻辑
        // （原代码中 cairo_eps_surface_get_file 是 EPS 专属函数，Cairo 1.18 已删除，必须移除）
        
        // 3. 抛出错误（明确提示 PDF 格式的上下文创建失败）
        throw std::runtime_error("saveRadialPowerSpectrum: failed to create cairo context for " + ext);
    }

        // -------------------------- 矢量图绘图配置（适配PDF） --------------------------
        // 1. 坐标变换：应用边距，确保内容不贴边
        cairo_translate(cr, margin, margin); // 平移到边距起点
        const Float content_w = plot_w - 2 * margin; // 有效绘图宽度（去除边距）
        const Float content_h = plot_h - 2 * margin; // 有效绘图高度
        cairo_scale(cr, 1.0, 1.0); // 矢量图无需缩放，保持原始尺寸精度

        // 2. 设置功率谱曲线样式（红色，线宽2pt，抗锯齿）
        cairo_set_source_rgb(cr, 1.0, 0.0, 0.0); // 红色（RGB）
        cairo_set_line_width(cr, 2.0); // 线宽2pt（矢量图单位，缩放后仍清晰）
        cairo_set_antialias(cr, CAIRO_ANTIALIAS_BEST); // 最佳抗锯齿，提升曲线平滑度

        // 3. 绘制功率谱曲线（基于计算的坐标点）
        for (int bin_idx = 0; bin_idx < num_bins - 1; bin_idx++) {
            // 计算矢量图坐标（x轴：径向距离区间，y轴：功率密度，适配有效绘图区）
            const Float x1 = static_cast<Float>(bin_idx) * bin_width;
            const Float x2 = static_cast<Float>(bin_idx + 1) * bin_width;
            // y轴翻转：Cairo 原点在左上角，需将"功率高"的点映射到上方
            const Float y1 = content_h - spectrum[bin_idx] * y_scale;
            const Float y2 = content_h - spectrum[bin_idx + 1] * y_scale;

            // 绘制线段（首次点用 move_to，后续用 line_to 连接）
            if (bin_idx == 0) {
                cairo_move_to(cr, x1, y1); // 起始点：移动画笔
            }
            cairo_line_to(cr, x2, y2); // 连接到下一个点
        }
        cairo_stroke(cr); // 执行曲线绘制（生效红色线条）

        // 4. 绘制坐标轴与网格（黑色线条，线宽1pt，虚线网格）
        cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); // 坐标轴颜色：黑色
        cairo_set_line_width(cr, 1.0); // 坐标轴线宽1pt

        // 4.1 绘制x轴（下边界）与y轴（左边界）
        cairo_move_to(cr, 0, content_h); // x轴起点（左下）
        cairo_line_to(cr, content_w, content_h); // x轴终点（右下）
        cairo_move_to(cr, 0, 0); // y轴起点（左上）
        cairo_line_to(cr, 0, content_h); // y轴终点（左下）
        cairo_stroke(cr); // 执行坐标轴绘制

        // 4.2 绘制网格线（灰色虚线，线宽0.5pt）
        cairo_set_source_rgb(cr, 0.8, 0.8, 0.8); // 网格颜色：浅灰色
        cairo_set_line_width(cr, 0.5); // 网格线宽0.5pt（不干扰主曲线）
        // 设置虚线样式：实线10pt + 空白5pt（适配矢量图缩放）
        const double dash_pattern[] = {10.0, 5.0};
        cairo_set_dash(cr, dash_pattern, 2, 0.0);

        // 绘制垂直网格线（每5个区间一条）
        const int grid_step = 5;
        for (int bin_idx = 0; bin_idx < num_bins; bin_idx += grid_step) {
            const Float x = static_cast<Float>(bin_idx) * bin_width;
            cairo_move_to(cr, x, 0); // 网格线起点（上边界）
            cairo_line_to(cr, x, content_h); // 网格线终点（下边界）
        }
        cairo_stroke(cr); // 执行网格绘制
        cairo_set_dash(cr, nullptr, 0, 0.0); // 重置为实线样式

        // 5. 添加坐标轴标签（矢量文本，适配边距位置）
        cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); // 文本颜色：黑色
        cairo_select_font_face(cr, "Arial", CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_NORMAL);
        cairo_set_font_size(cr, 12.0); // 字体大小12pt（矢量文本，缩放后清晰）

        // 5.1 x轴标签："Normalized Radial Distance"（位于x轴下方，居中）
        cairo_text_extents_t x_text_extents;
        const std::string x_label = "Normalized Radial Distance";
        cairo_text_extents(cr, x_label.c_str(), &x_text_extents);
        const Float x_label_x = (content_w - x_text_extents.width) / 2.0; // 水平居中
        const Float x_label_y = content_h + margin - 10.0; // 位于x轴下方10pt
        cairo_move_to(cr, x_label_x, x_label_y);
        cairo_show_text(cr, x_label.c_str());

        // 5.2 y轴标签："Power Density"（垂直旋转，位于y轴左侧，居中）
        const std::string y_label = "Power Density";
        cairo_save(cr); // 保存当前绘图状态（避免旋转影响后续内容）
        // 旋转文本：绕y轴中点逆时针旋转90度
        const Float y_label_rotate_x = -10.0; // 位于y轴左侧10pt
        const Float y_label_rotate_y = content_h / 2.0; // 垂直居中
        cairo_translate(cr, y_label_rotate_x, y_label_rotate_y);
        cairo_rotate(cr, -M_PI / 2.0); // 逆时针旋转90度（弧度制）
        // 计算垂直文本的居中位置
        cairo_text_extents(cr, y_label.c_str(), &x_text_extents);
        const Float y_label_text_y = (content_h + x_text_extents.width) / 2.0;
        cairo_move_to(cr, 0, y_label_text_y);
        cairo_show_text(cr, y_label.c_str());
        cairo_restore(cr); // 恢复绘图状态（取消旋转）

        // -------------------------- 释放 Cairo 资源（避免内存泄漏） --------------------------
        cairo_show_page(cr); // 刷新页面（确保矢量图内容写入文件）
        cairo_destroy(cr); // 销毁绘图上下文
        cairo_surface_finish(surface); // 完成表面写入（强制刷写缓冲区）
        cairo_surface_destroy(surface);

    } else {
        // 不支持的文件格式：抛出异常并提示支持的格式
        throw std::invalid_argument(
            "saveRadialPowerSpectrum: unsupported format - " + ext + 
            " (only .png, .pdf are supported)"
        );
    }

    std::cout << "[Info] Radial power spectrum saved to: " << file_path << std::endl;
}


// 关键：若实现放在 .cpp 文件中，必须显式实例化 float 版本（否则链接器找不到）
template class GBNBasePointSet<float>; // 显式实例化 float 模板，生成对应成员函数代码